/*
 * Example_CalciteMagnesiteDolomite.h
 *
 *  Created on: 19 Oct 2011
 *      Author: allan
 */

// C++ includes
#include <fstream>
#include <iostream>
#include <iomanip>
using namespace std;

// GeoReact includes
#include <GeoReact.h>

//----------------------------------------------------------------------//
// Rock Composition:
// 
//     65% of Calcite   (100.09 g/mol, volumic mass: 2709.89 kg/m3)
//     25% of Magnesite ( 84.31 g/mol, volumic mass: 3009.28 kg/m3)
//     10% of Dolomite  (184.40 g/mol, volumic mass: 2866.04 kg/m3)
// 
// Total mass of rock: 200 g
// 
// Mass of Calcite:   130 g
// Mass of Magnesite:  50 g
// Mass of Dolomite:   20 g
// 
// Moles of Calcite:   1.29883105 mol
// Moles of Magnesite: 0.59304946 mol
// Moles of Dolomite:  0.10845987 mol
//----------------------------------------------------------------------//

void Example_CalciteMagnesiteDolomite()
{
	const double nCalcite   = 130.0 / 100.09; // conversion from g to mol
	const double nMagnesite =  50.0 /  84.31; // conversion from g to mol
	const double nDolomite  =  20.0 / 184.40; // conversion from g to mol

	const double ssaCalcite   = 50.0 * 100.09 * 1.0E-4; // conversion from cm2/g to m2/mol
	const double ssaMagnesite = 50.0 *  84.31 * 1.0E-4; // conversion from cm2/g to m2/mol
	const double ssaDolomite  = 50.0 * 184.40 * 1.0E-4; // conversion from cm2/g to m2/mol

	GeoManager manager;
	
	manager.SetAqueousSpecies({"H2O", "Ca[2+]", "HCO3[-]", "H[+]", "Mg[2+]", "CO2(aq)", "CO3[2-]", "CaCO3(aq)", "CaHCO3[+]", "CaOH[+]", "MgCO3(aq)", "MgHCO3[+]", "OH[-]"});//, "Mg4(OH)4[4+]"});
	manager.SetMineralSpecies({"Calcite", "Magnesite", "Dolomite"});
	
	MineralKineticsInfo kCalcite("Calcite");
	kCalcite.AddMineralMechanism(pow(10, -5.81), 23.5);
	kCalcite.SetSpecificSurfaceArea(ssaCalcite);
	
	MineralKineticsInfo kMagnesite("Magnesite");
	kMagnesite.AddMineralMechanism(pow(10, -9.34), 23.5);
	kMagnesite.SetSpecificSurfaceArea(ssaMagnesite);
	
	MineralKineticsInfo kDolomite("Dolomite");
	kDolomite.AddMineralMechanism(pow(10, -7.53), 52.2);
	kDolomite.SetSpecificSurfaceArea(ssaDolomite);
	
	manager.SetMineralKinetics("Calcite", kCalcite);
	manager.SetMineralKinetics("Magnesite", kMagnesite);
	manager.SetMineralKinetics("Dolomite", kDolomite);
	
	GeoSystem system = manager.AssembleGeoSystem();
	
	cout << system << endl;
	
	GeoEquilibrator equilibrator(system);
	
	equilibrator.ImposeChargeBalance();
	equilibrator.ImposeSpeciesMass("H2O", 55.508);
	equilibrator.ImposeSpeciesMass("Ca[2+]", 1.0E-7);
	equilibrator.ImposeSpeciesMass("Mg[2+]", 1.0E-7);
	equilibrator.ImposeSpeciesMass("HCO3[-]", 1.0E-7);
	equilibrator.ImposeSpeciesMass("Calcite", nCalcite);
	equilibrator.ImposeSpeciesMass("Magnesite", nMagnesite);
	equilibrator.ImposeSpeciesMass("Dolomite", nDolomite);
	
	GeoState state(system);
	
	equilibrator.Equilibrate(state);
	
	state.SetComposition(state.GetComposition().array().abs());
	
	cout << state << endl;
	
	GeoReactor reactor(system);
	
	reactor.SetAlgorithm(RungeKuttaFehlberg45);
	reactor.SetAdaptiveTimeStepControl(1.0E-5, 0);
	
	const double hour = 3600.0;
	const double day  = 24 * hour;
	const double year = 365 * day;
	
	double t = 0.0, tEnd = 10 * year;
	
	ofstream file("result.dat");
	
	cout << setprecision(10);
	file << setprecision(10);
	
	reactor.SetMinimumTimeStep(day);
	
	while(t < tEnd)
	{
		cout << setw(20) << left << t/day;
		cout << setw(20) << left << state.Get("Moles", "Calcite") * 100.09;
		cout << setw(20) << left << state.Get("Moles", "Magnesite") * 84.31;
		cout << setw(20) << left << state.Get("Moles", "Dolomite") * 184.40;
		cout << setw(20) << left << state.Get("Moles", "Ca[2+]");
		cout << setw(20) << left << state.Get("Moles", "Mg[2+]");
		cout << setw(20) << left << state.Get("Moles", "HCO3[-]");
		cout << setw(20) << left << state.Get("Moles", "CO2(aq)");
		cout << setw(20) << left << state.Get("KineticRate", "Calcite");
		cout << setw(20) << left << state.Get("KineticRate", "Ca[2+]");
		cout << setw(20) << left << state.Get("pH") << endl;
		
		file << setw(20) << left << t/day;
		file << setw(20) << left << state.Get("Moles", "Calcite") * 100.09;
		file << setw(20) << left << state.Get("Moles", "Magnesite") * 84.31;
		file << setw(20) << left << state.Get("Moles", "Dolomite") * 184.40;
		file << setw(20) << left << state.Get("Moles", "Ca[2+]");
		file << setw(20) << left << state.Get("Moles", "Mg[2+]");
		file << setw(20) << left << state.Get("Moles", "HCO3[-]");
		file << setw(20) << left << state.Get("Moles", "CO2(aq)");
		file << setw(20) << left << state.Get("KineticRate", "Calcite");
		file << setw(20) << left << state.Get("KineticRate", "Ca[2+]");
		file << setw(20) << left << state.Get("pH") << endl;
		
		reactor.ReactEuler(t, tEnd, state);
		//reactor.React(t, tEnd, state);
	}
}

void Example_CalciteMagnesiteDolomiteEquilibrium()
{
	const double nCalcite   = 130.0 / 100.09; // conversion from g to mol
	const double nMagnesite =  50.0 /  84.31; // conversion from g to mol
	const double nDolomite  =  20.0 / 184.40; // conversion from g to mol

	GeoManager manager;
	
	manager.SetAqueousSpecies({"H2O", "Ca[2+]", "HCO3[-]", "H[+]", "Mg[2+]", "CO2(aq)", "CO3[2-]", "CaCO3(aq)", "CaHCO3[+]", "CaOH[+]", "MgCO3(aq)", "MgHCO3[+]", "OH[-]"});//, "Mg4(OH)4[4+]"});
	//manager.SetGaseousSpecies({"H2O(g)", "CO2(g)"});
	manager.SetMineralSpecies({"Calcite", "Magnesite", "Dolomite"});
	
	GeoSystem system = manager.AssembleGeoSystem();
	
	cout << system << endl;
	
	GeoEquilibrator equilibrator(system);
	
	equilibrator.ImposeChargeBalance();
	equilibrator.ImposeSpeciesMass("H2O", 55.508);
	equilibrator.ImposeElementMass("Ca", nCalcite + nDolomite);
	equilibrator.ImposeElementMass("Mg", nMagnesite + nDolomite);
	equilibrator.ImposeElementMass("C", nCalcite + nMagnesite + 2 * nDolomite);
	
	GeoState state(system);
	
	equilibrator.Equilibrate(state);
	
	cout << state << endl;
}



